****************************************************************************
*``The Costs and Environmental Justice Concerns of NIMBY in Solid Waste Disposal''*
* Phuong Ho (SNF-Centre for Applied Research in Norwegian School of Economics) *
*******************************************************************************
*** This file makes Figure 1, Figure 2, Table 1, Table 2 (first stage results),
** Table B1, Figure B1, Figure C1, Figure C2
*The file also generates inputs for the structural models in nimby_main.m
cls
clear all 
capture log close
set more off

global mydir "C:\XYZ" // set working directory
global data "$mydir\data"
global matlabin "$mydir\matlab\input"
global result "$mydir\TEX\"
cd "$mydir\stata"
set matsize 11000
set maxvar 32700

*************** IMPORT RAW MASTER DATA ***************************************
use "$data/master_data.dta", clear
**This data set is the raw data, i.e. has not removed any observations. 
**The data just combined all trash flow data and matched code from...
*...Waste Business Journal (WBJ) data for disposal fee (price) information.
*Note: the data also include pairs of county-facility that have zero amount of trash
*which may include facilities, at a time t, that have total trash over counties=0
******************************************************************************
** NOW, CLEAN DATA FOR ANALYSIS / REMOVE SHUT-DOWN AND UNMATCHED OBSERVATIONS **
***Remove those shut-down facilities***
bys year quarter: egen CAq0=total(q) // total trash generated by CA
sort year quarter county facility
by year quarter county: egen cntyq0=total(q) // total trash generated by a county
sort year quarter facility
by year quarter facility: egen totq_f = total(q) // total trash at a facility
tab nocode [iw=q] // 0.52% of trash at facilities not in price dataset
tab zerop [iw=q] // 0.41% of trash at facilities with price=0
tab island [iw=q] // 0.01% of trash at facilities in islands (in LA county)
tab exp [iw=q] // 1.16% of trash exported to outside CA
tab county exp [iw=q], row nofreq 
* remove those lack-of-information observations
keep if totq_f >0 & nocode==0 & zerop==0 &island==0

*for those that have the same code, keep the one with shortest distance
sort year quarter county Code facility
by year quarter county Code: egen quan = total(q)
replace q = quan 
drop quan
sort year quarter county Code dridis 
gen mindist = 0
by year quarter county Code: replace mindist = 1 if _n==1
keep if mindist==1

gen j= Code
replace j = "oos" if facility=="oos"
replace descnty = "Yuba-Sutter" if descnty=="Yuba"
sort year
merge m:1 year using index2000
drop if _merge==2
drop _merge
gen realprice = price * index2000
label var realprice "Price in 2000 dollars"

sort year quarter county j

*************SAVE TO DATA FOR ANALYSIS ********************************
keep if year <= 2015
keep year quarter county facility exp q Code j price  ///
	posq expcnty cpwlat cpwlon deslat deslon descnty ///
	oos Vdis dridis yq index2000 realprice cntyq0 CAq0 ///
	Operation OperatorEntity OperatorCode Ownership OwnerEntity OwnerCode ///
	Access DaysYear cntyland cntyintlat cntyintlon
gen localf = county==descnty
gen radius = .

***** generate categorical variables of groups *****
sort county
egen C = group(county), label
gen desco=descnty
replace desco="oos" if facility=="oos"
sort desco
egen D = group(desco), label
sort county desco
egen CD=group(C D), label
sort Code
egen F = group(Code), label
sort j
egen J = group(j), label
sort county j
egen CJ = group(county j)
sort year quarter
egen T = group(year quarter)
sort T county
egen TC = group(T county)
sort T j
egen TJ = group(T j)
sort T Code
egen TF = group(T Code), label
sort T desco
egen TD=group(T desco)
sort TC j
egen TCJ = group(TC j)

save nimby_main, replace // 673,968 obs

******************************************************************************
*************** GRAPHICAL DESCRIPTIVE STATS ********************************
**************************************************************************
**************** Figure 1F ************************************************
use nimby_main, clear
keep if q>0
local rad 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150
gen rad = .
gen outlierqpc = .
tokenize `rad'
local i 1
while "``i''" != "" {
	qui replace rad = ``i'' in `i'
	replace radius = ``i''
	sort yq county Code
    gen outlier = dridis > radius
	egen sumq = total(q)
	egen outliq_x = total(q) if outlier==1
	egen outliq = mean(outliq_x)
	replace outlierqpc = outliq/sumq*100 in `i'
	drop outlier sumq outliq_x outliq
	local ++i
}
sum dridis
local r = r(max)
local rr = round(`r',.01)
twoway line  outlierqpc rad, ///
	   	lpattern("solid") lwidth(medthick) lcolor(black) ///
		legend(off) ///
		graphregion(color(white)) ///
		ylabel(#10, grid) ///
		xmtick(20 30 40 60 70 80 90 110 120 130 140 150) ///
		xtick() ///
		xlabel(10 30 40 50 60 100 150, grid) c(1) yaxis(1) ///
		xtitle("miles",placement("east")) ytitle("% trash disposed of outside a given mile", placement("north"))
graph export "$result/figure1f.pdf", replace

******************************* FIGURE 1 *************************************
use nimby_main, clear
sort TC
by TC: egen tc_q = total(q)
by TC: egen inq_del=total(q) if localf==1
by TC: egen inq=mean(inq_del)
gen inq_p=inq*100/tc_q
by TC: egen ouq_del=total(q) if localf==0
by TC: egen ouq=mean(ouq_del)
gen ouq_p=ouq*100/tc_q
by TC: gen tcdup=cond(_N==1,0,_n)

local varl "tc_q inq_p ouq_p inq ouq"
sort T
foreach v in `varl'{
	by T: egen `v'_t_del = mean(`v') if tcdup<=1
	by T: egen `v'_t=mean(`v'_t)
}
drop *_del
label var inq_p_t "Percentage of waste a county dispose locally"
label var ouq_p_t "Percentage of waste a county sends to other counties"
label var inq_t "Amount of waste a county dispose locally"
label var ouq_t "Amount of waste a county sends to other counties"
label var tc_q_t "Amount of waste a county generates"

replace inq_t=inq_t/1000
replace ouq_t=ouq_t/1000
replace tc_q_t=tc_q_t/1000

by T: egen t_q=total(q)

gen disw=q*dridis
by T: egen t_disw=total(disw)
gen t_dis=t_disw/t_q
label var t_dis "(Waste-weighted) Distance a county transports waste"

gen priw=q*realprice
by T: egen t_priw=total(priw)
gen t_pri=t_priw/t_q
label var t_pri "(Waste-weighted) disposal price (in 2000 \$) haulers pay"

**number of facilities
keep if  facility!="oos" // only count # facilities within California
bys yq Code: gen dup= cond(_N==1,0,_n)
keep if dup<=1
sort T
by T: egen noj = count(j)

by T: gen tdup=cond(_N==1,0,_n)
keep if tdup<=1
keep if year<=2015
************************************ Figure 1a ********************************
twoway (line tc_q_t yq, ///
       lpattern("solid" ) lwidth(medthick) lcolor(navy black) ///
	   	   legend(off) ///
		graphregion(color(white)) ///
		ylabel(#10, grid) ///
		xlabel(#10, grid angle(45)) c(1) yaxis(1) ///
		xtitle("",placement("east")) ytitle("thousand tons", placement("north") axis(1)) ) 
graph export "$result/figure1a.pdf", replace

***************************** Figure 1b **************************************
twoway (line noj yq, ///
       lpattern("solid" ) lwidth(medthick) lcolor(navy black) ///
	   	   legend(off) ///
		graphregion(color(white)) ///
		ylabel(#10, grid) ///
		xlabel(#10, grid angle(45)) c(1) yaxis(1) ///
		xtitle("",placement("east")) ytitle("# facilities", placement("north") axis(1)) ) 
graph export "$result/figure1b.pdf", replace

********************************** Figure 1c *********************************
twoway (line t_pri yq, ///
       lpattern("solid" ) lwidth(medthick) lcolor(navy black) ///
	   	   legend(off) ///
		graphregion(color(white)) ///
		ylabel(#10, grid) ///
		xlabel(#10, grid angle(45)) c(1) yaxis(1) ///
		xtitle("",placement("east")) ytitle("$/ton", placement("north") axis(1)) ) 
graph export "$result/figure1c.pdf", replace

********************************* Figure 1d ***********************************
twoway (line ouq_p_t yq, ///
       lpattern("solid" ) lwidth(medthick) lcolor(navy black) ///
	   	   legend(off) ///
		graphregion(color(white)) ///
		ylabel(#10, grid) ///
		xlabel(#10, grid angle(45)) c(1) yaxis(1) ///
		xtitle("",placement("east")) ytitle("export as a percentage of a county waste (%)", placement("north") axis(1)) ) 
graph export "$result/figure1d.pdf", replace

********************************* Figure 1e ***********************************
twoway (line t_dis yq, ///
       lpattern("solid" ) lwidth(medthick) lcolor(navy black) ///
	   	   legend(off) ///
		graphregion(color(white)) ///
		ylabel(#10, grid) ///
		xlabel(#10, grid angle(45)) c(1) yaxis(1) ///
		xtitle("",placement("east")) ytitle("mile", placement("north") axis(1)) ) 
graph export "$result/figure1e.pdf", replace

*****************************************************************************
*DEFINING CHOICE SET. EXAMINE PRICE RESPONSE AND DISTANCE RESPONSE BY DISTANCE
******##############################################################*********
*********************** TABLE B1 AND FIGURE B1 *******************************
****************************************************************************
************************* Trash flows within 120 miles *********************
************* Columns (1), (2) in table B1. And figures B1-A, B1-C ***********
use nimby_main, clear
replace radius = 120 // set the radius that defines a choice set
sort yq county Code
gen outlier = dridis > radius
**Calculate market share of facility in each market
cap gen mktshare00=q/cntyq0*100 if outlier==0

** generate linear splines for dridis
keep if dridis<=radius
mkspline dridis1 30 dridis2 60 dridis3 90 dridis4 = dridis 
local ni 4
local knots "0 30 60 90"
tokenize `knots'
local i 1
forval i=1/`ni' {
	if (`i' < `ni') {
		local k = `i'+1
		gen dridisk`i' = (dridis >=``i'' & dridis < ``k'')
	}
	if (`i' == `ni') {
		gen dridisk`i' = (dridis >=``i'')
	}
} 
local varl "realprice"
foreach v in `varl' {
  forval i=1/`ni' {
	gen `v'`i' = `v'*dridisk`i'
	*local ++i
  }
}

label var realprice "price"
label var dridis "distance"

sum radius
local r = r(max)
local rr = round(`r',1)

eststo clear
lsfereg mktshare00 realprice1-realprice`ni' dridis1-dridis`ni', cluster(C) absorb(TC J)
estimates store mod1
lsfereg mktshare00 realprice1-realprice`ni', cluster(C) absorb(T CJ)
estimates store mod2
label var dridis1 "distance:0-30"
label var dridis2 "distance:30-60"
label var dridis3 "distance:60-90"
label var dridis4 "distance:90-120"
label var realprice1 "price:0-30"
label var realprice2 "price:30-60"
label var realprice3 "price:60-90"
label var realprice4 "price:90-120"
***********  COLUMNS (1) AND (2) OF TABLE B1 *************************
esttab mod1 mod2 using "$result/tableB1_12.tex" , replace ///
	label  se(4) ar2(4) nonotes scalars(N_unc ll) compress b(4) nogaps  ///
		title("Distance response and price response by distance, r = `rr'") ///
		star(* 0.10 ** 0.05 *** 0.01)  ///
		keep(*price* dridis*) ///
		prefoot(\hline "facility FE &  Y & Y "  \\ ///
		"quarter FE & Y  & Y " \\ ///
		"quarter x origin cnty FE & Y  &  " \\ ///
		"origin cnty x facility FE &   & Y " \\ ///
		"Cluster SE & origin cnty & origin cnty " \\ \hline ) 

label var dridis1 "distance:0-30"
label var dridis2 "30-60"
label var dridis3 "60-90"
label var dridis4 "90-120"
label var realprice1 "distance:0-30"
label var realprice2 "30-60"
label var realprice3 "60-90"
label var realprice4 "90-120"
************************** FIGURE B1-A **************************************
coefplot mod1, drop(_cons *price*) yline(0) vertical xsize(4) graphregion(color(none)) bgcolor(white)
graph export "$result/figureB1a.pdf", replace
************************** FIGURE B1-C **************************************
coefplot mod2, drop(_cons) yline(0) vertical xsize(4) graphregion(color(none)) bgcolor(white)
graph export "$result/figureB1c.pdf", replace

************************** Trash flows within 150 miles *********************
******************* Columns (3), (4) in table B1 *****************************
use nimby_main, clear
replace radius = 150 // set the radius that defines a choice set
sort yq county Code
gen outlier = dridis > radius
**Calculate market share of facility in each market
cap gen mktshare00=q/cntyq0*100 if outlier==0

** generate linear splines for dridis
keep if dridis<=radius
mkspline dridis1 40 dridis2 70 dridis3 100 dridis4 130 dridis5 = dridis 
local ni 5
local knots "0 40 70 100 130"
tokenize `knots'
local i 1
forval i=1/`ni' {
	if (`i' < `ni') {
		local k = `i'+1
		gen dridisk`i' = (dridis >=``i'' & dridis < ``k'')
	}
	if (`i' == `ni') {
		gen dridisk`i' = (dridis >=``i'')
	}
} 
local varl "realprice"
foreach v in `varl' {
  forval i=1/`ni' {
	gen `v'`i' = `v'*dridisk`i'
	*local ++i
  }
}

label var realprice "price"
label var dridis "distance"

sum radius
local r = r(max)
local rr = round(`r',1)

eststo clear
lsfereg mktshare00 realprice1-realprice`ni' dridis1-dridis`ni', cluster(C) absorb(TC J)
estimates store mod1
lsfereg mktshare00 realprice1-realprice`ni', cluster(C) absorb(T CJ)
estimates store mod2
label var dridis1 "distance:0-40"
label var dridis2 "distance:40-70"
label var dridis3 "distance:70-100"
label var dridis4 "distance:100-130"
label var dridis5 "distance:130-150"
label var realprice1 "price:0-40"
label var realprice2 "price:40-70"
label var realprice3 "price:70-100"
label var realprice4 "price:100-130"
label var realprice5 "price:130-150"
****************  COLUMNS (3) AND (4) OF TABLE B1 *************************
esttab mod1 mod2 using "$result/tableB1_34.tex" , replace ///
	label  se(4) ar2(4) nonotes scalars(N_unc ll) compress b(4) nogaps  ///
		title("Distance response and price response by distance, r = `rr'") ///
		star(* 0.10 ** 0.05 *** 0.01)  ///
		keep(*price* dridis*) ///
		prefoot(\hline "facility FE &  Y & Y "  \\ ///
		"quarter FE & Y  & Y " \\ ///
		"quarter x origin cnty FE & Y  &  " \\ ///
		"origin cnty x facility FE &   & Y " \\ ///
		"Cluster SE & origin cnty & origin cnty " \\ \hline ) 

************************** Trash flows within 150 miles *********************
********** Columns (5), (6) in table B1. And figures B1-B, B1-D *************
use nimby_main, clear
replace radius = 150 // set the radius that defines a choice set
sort yq county Code
gen outlier = dridis > radius
**Calculate market share of facility in each market
cap gen mktshare00=q/cntyq0*100 if outlier==0

** generate linear splines for dridis
keep if dridis<=radius
mkspline dridis1 30 dridis2 60 dridis3 90 dridis4 110 dridis5 130 dridis6 = dridis 
local ni 6
local knots "0 30 60 90 110 130"
tokenize `knots'
local i 1
forval i=1/`ni' {
	if (`i' < `ni') {
		local k = `i'+1
		gen dridisk`i' = (dridis >=``i'' & dridis < ``k'')
	}
	if (`i' == `ni') {
		gen dridisk`i' = (dridis >=``i'')
	}
} 
local varl "realprice"
foreach v in `varl' {
  forval i=1/`ni' {
	gen `v'`i' = `v'*dridisk`i'
	*local ++i
  }
}

label var realprice "price"
label var dridis "distance"

sum radius
local r = r(max)
local rr = round(`r',1)

eststo clear
lsfereg mktshare00 realprice1-realprice`ni' dridis1-dridis`ni', cluster(C) absorb(TC J)
estimates store mod1
lsfereg mktshare00 realprice1-realprice`ni', cluster(C) absorb(T CJ)
estimates store mod2
label var dridis1 "distance:0-30"
label var dridis2 "distance:30-60"
label var dridis3 "distance:60-90"
label var dridis4 "distance:90-110"
label var dridis5 "distance:110-130"
label var dridis6 "distance:130-150"
label var realprice1 "price:0-30"
label var realprice2 "price:30-60"
label var realprice3 "price:60-90"
label var realprice4 "price:90-110"
label var realprice5 "price:110-130"
label var realprice6 "price:130-150"
****************  COLUMNS (5) AND (6) OF TABLE B1 *************************
esttab mod1 mod2 using "$result/tableB1_56.tex" , replace ///
	label  se(4) ar2(4) nonotes scalars(N_unc ll) compress b(4) nogaps  ///
		title("Distance response and price response by distance, r = `rr'") ///
		star(* 0.10 ** 0.05 *** 0.01)  ///
		keep(*price* dridis*) ///
		prefoot(\hline "facility FE &  Y & Y "  \\ ///
		"quarter FE & Y  & Y " \\ ///
		"quarter x origin cnty FE & Y  &  " \\ ///
		"origin cnty x facility FE &   & Y " \\ ///
		"Cluster SE & origin cnty & origin cnty " \\ \hline ) 

label var dridis1 "dist.:0-30"
label var dridis2 "30-60"
label var dridis3 "60-90"
label var dridis4 "90-110"
label var dridis5 "110-130"
label var dridis6 "130-150"
label var realprice1 "dist.:0-30"
label var realprice2 "30-60"
label var realprice3 "60-90"
label var realprice4 "90-110"
label var realprice5 "110-130"
label var realprice6 "130-150"
************************** FIGURE B1-B **************************************
coefplot mod1, drop(_cons *price*) yline(0) vertical xsize(4) graphregion(color(none)) bgcolor(white)
graph export "$result/figureB1b.pdf", replace
************************** FIGURE B1-D **************************************
coefplot mod2, drop(_cons) yline(0) vertical xsize(4) graphregion(color(none)) bgcolor(white)
graph export "$result/figureB1d.pdf", replace

************************** Trash flows within 150 miles *********************
********** Columns (7), (8) in table B1 ***************
use nimby_main, clear
replace radius = 150 // set the radius that defines a choice set
sort yq county Code
gen outlier = dridis > radius
**Calculate market share of facility in each market
cap gen mktshare00=q/cntyq0*100 if outlier==0

** generate linear splines for dridis
keep if dridis<=radius
local r = r(max)
local rr = round(`r',1)
mkspline dridis1 30 dridis2 60 dridis3 80 dridis4 100 dridis5 125 dridis6 = dridis 
local ni 6
local knots "0 30 60 80 100 125"
tokenize `knots'
local i 1
forval i=1/`ni' {
	if (`i' < `ni') {
		local k = `i'+1
		gen dridisk`i' = (dridis >=``i'' & dridis < ``k'')
	}
	if (`i' == `ni') {
		gen dridisk`i' = (dridis >=``i'')
	}
} 
local varl "realprice"
foreach v in `varl' {
  forval i=1/`ni' {
	gen `v'`i' = `v'*dridisk`i'
	*local ++i
  }
}

label var realprice "price"
label var dridis "distance"

sum radius
local r = r(max)
local rr = round(`r',1)

eststo clear
lsfereg mktshare00 realprice1-realprice`ni' dridis1-dridis`ni', cluster(C) absorb(TC J)
estimates store mod1
lsfereg mktshare00 realprice1-realprice`ni', cluster(C) absorb(T CJ)
estimates store mod2
label var dridis1 "distance:0-30"
label var dridis2 "distance:30-60"
label var dridis3 "distance:60-80"
label var dridis4 "distance:80-100"
label var dridis5 "distance:100-125"
label var dridis6 "distance:125-150"
label var realprice1 "price:0-30"
label var realprice2 "price:30-60"
label var realprice3 "price:60-80"
label var realprice4 "price:80-100"
label var realprice5 "price:100-125"
label var realprice6 "price:125-150"
****************  COLUMNS (7) AND (8) OF TABLE B1 *************************
esttab mod1 mod2 using "$result/tableB1_78.tex" , replace ///
	label  se(4) ar2(4) nonotes scalars(N_unc ll) compress b(4) nogaps  ///
		title("Distance response and price response by distance, r = `rr'") ///
		star(* 0.10 ** 0.05 *** 0.01)  ///
		keep(*price* dridis*) ///
		prefoot(\hline "facility FE &  Y & Y "  \\ ///
		"quarter FE & Y  & Y " \\ ///
		"quarter x origin cnty FE & Y  &  " \\ ///
		"origin cnty x facility FE &   & Y " \\ ///
		"Cluster SE & origin cnty & origin cnty " \\ \hline ) 

*******************************************************************************
***************** FIGURE F2 IN APPENDIX F *********************************
********** Export into csv data to plot a map of county and facility **********
use nimby_main, clear
*total waste generated by a county
bys year county: egen cnty_generateq=total(q)
*total waste a county exports to other counties
bys year county: egen cnty_exportq_del=total(q) if county~=desco
bys year county: egen cnty_exportq=mean(cnty_exportq_del)
replace cnty_exportq=0 if cnty_exportq==.
drop *_del
*total waste a county receives
bys year desco: egen cnty_receiveq_del=total(q)
gen cnty_receiveq_del2=cnty_receiveq_del if county==desco // spread values
bys year county: egen cnty_receiveq=mean(cnty_receiveq_del2)
replace cnty_receiveq=0 if cnty_receiveq==.
drop *_del*
*total imported waste a county receives
bys year desco: egen cnty_importq_del =total(q) if county~=desco & j~="oos" 
bys year desco: egen cnty_importq_del2=mean(cnty_importq_del) // spread values to des cnty level
replace cnty_importq_del2=0 if cnty_importq_del2==.
gen cnty_importq_del3=cnty_importq_del2 if county==desco // spread values from des cnty to county
bys year county: egen cnty_importq=mean(cnty_importq_del3) // spread values to county
replace cnty_importq=0 if cnty_importq==.
drop *_del*
*total waste a facility receives. Consider oos is 1 facility, but we won't show it on map
bys year facility: egen fcty_q=total(q)
*total waste a facility receives from counties other than its home county
bys year facility: egen fcty_importq_del=total(q) if county~=desco
bys year facility: egen fcty_importq=mean(fcty_importq_del)
drop *_del
bys year county: gen dupc = cond(_N==1,0,_n)
bys year facility: gen dupj=cond(_N==1,0,_n)

preserve
keep if oos==0
bys facility: gen dup=cond(_N==1,0,_n)
keep if dup<=1
keep facility Code deslat deslon j
export delimited using "$result/map_facility_location_1995_2015", nolabel replace 
restore
preserve
keep if year==2010
keep if dupc<=1
keep year county cpwlat cpwlon cnty_generateq cnty_exportq cnty_receiveq cnty_importq
expand 2 if county=="Yuba-Sutter" // replace obs if by 2 copies of it
replace county="Yuba" if county=="Yuba-Sutter"
replace county="Sutter" if _n==_N // last observation
export delimited using "$result/map_county_waste", nolabel replace 
restore
preserve
keep if year==2010
keep if dupj<=1
keep year facility deslat deslon descnty price realprice Code fcty_q fcty_importq
export delimited using "$result/map_facility_waste_2010", nolabel replace 
restore

****************************************************************************
**** GENERATE DATA SETS TO ESTIMATE THE STRUCTURAL MODEL *********
**** PREPARE INSTRUMENTS AND FIRST STAGE RESULTS *********
use nimby_main, clear
replace radius = 60 
sort TC Code
gen outlier = dridis > radius
tab exp [iw=q] if outlier==0
gen share=q/cntyq0
gen share00=share*100
********** Calculate the market size of each county ******************
** Relevant market size in each area/county = sum over facilities in the market
sort TC Code
by TC: egen mktsize_del = total(q) if outlier==0
by TC: egen mktsize = mean(mktsize_del) // to spread the values
drop mktsize_del
gen outq0=cntyq0-mktsize
sort yq Code
by yq Code: egen summktsizes_del = total(mktsize) if outlier==0
by yq Code: egen summktsizes = mean(summktsizes_del) 

* Total market sizes of other competing markets of a facility at t
* = summktsizes - mktsize of a specific competing market
gen sumothmktsizes = summktsizes - mktsize 
gen sumothmktsizes50s = sumothmktsizes/100000 

*** Number of competing markets for a facility at a time t
sort yq Code
by yq Code: egen nummkt_del = count(county) if outlier==0 
by yq Code: egen nummkt = mean(nummkt_del) 
gen numothmkt = nummkt - 1 
gen isumothmktsizes50s = (numothmkt > 0 & numothmkt < .)
replace isumothmktsizes50s = . if mktsize<0.01 // county that has all trash flows >r
drop *_del

***********##############################################################****
************ SUMMARY STATISTICS OF MAIN VARIABLES **************************
************************* TABLE 1 *******************************************
****************************** PANEL A2 AND B2 ******************************
preserve // summary statistics of ALL trash flows conditional on positive q
keep if q>0 
eststo clear
// Panel A2
estpost sum q dridis realprice
eststo mod1a
estpost sum dridis realprice [iw=q] // waste weighted distance and price
eststo mod1b

// Panel B2
sort yq county // unit: quarter x origin county
egen q_tc = total(q), by(yq county)
egen outq_tc_x = total(q) if localf==0, by(yq county)
egen outq_tc  = mean(outq_tc_x), by(yq county) // total trash goes outside county
replace outq_tc = 0 if outq_tc ==.
gen outqpc_tc = outq_tc*100/q_tc  // pecentage of trash goes outside
drop *_x
by yq county: gen dup = cond(_N==1,0,_n)
estpost sum q_tc outqpc_tc if dup <=1
eststo mod2
esttab mod1a mod1b mod2 using "$result/table1_a2b2.tex", cells("count mean(fmt(%9.2f)) sd(fmt(%9.2f)) min(fmt(%9.2f)) max") noobs replace
restore

****************************** PANEL A1 AND B1 ********************************
preserve // summary statistics of the balance sample for analysis. R = 60 miles
keep if outlier==0 & mktsize>0 
eststo clear
// Panel A1
estpost sum q dridis realprice
eststo mod1a
estpost sum dridis realprice [iw=q] // waste weighted distance and price
eststo mod1b

// Panel B1
sort yq county // unit: quarter x origin county
egen q_tc = total(q), by(yq county)
egen outq_tc_x = total(q) if localf==0, by(yq county)
egen outq_tc  = mean(outq_tc_x), by(yq county) // total trash goes outside county
replace outq_tc = 0 if outq_tc ==.
gen outqpc_tc = outq_tc*100/q_tc  // pecentage of trash goes outside
drop *_x
sort yq county j
egen tcj = group(yq county j)
sort yq county
egen nof_ct = count(tcj), by(yq county) // number of options in a market
by yq county: gen dup = cond(_N==1,0,_n)

estpost sum q_tc outqpc_tc nof_ct if dup <=1
eststo mod2
drop dup

esttab mod1a mod1b mod2 using "$result/table1_a1b1.tex", cells("count mean(fmt(%9.2f)) sd(fmt(%9.2f)) min(fmt(%9.2f)) max") noobs replace
restore

*****************************************************************************
**** GENERATE DATA SETS TO ESTIMATE THE STRUCTURAL MODEL *******************
keep if outlier==0 & mktsize>0 
sort TC Code
save nimby_estimate, replace
**************************************************************************
****** IMPORT DEMOGRAPHICS IN 3-MILE COMMUNITIES *************************
// 1/3. Import demographics for county of origin 
use $data\demographics\raw\county_demographics_2010_CA, clear
keep *geoid *name *pop *white* *black* *hispanic* *asian* *mdhhinc
foreach x of varlist _all {
	rename `x' O_`x'
} 
rename O_COname county
gen year=2010
sort year county
tempfile tmp
save `tmp', replace
use nimby_estimate, clear
sort year county
merge m:1 year county using `tmp'
erase `tmp'
drop if _merge==2 //_merge==1 are not 2010, _merge==2 are counties that doesn't generate trash in 2010
drop _merge
replace O_COmdhhinc=O_COmdhhinc*index2000
save nimby_estimate, replace

// 2/3. Import demographics for county of destination 
use $data\demographics\raw\county_demographics_2010_CA, clear
keep *geoid *name *pop *white* *black* *hispanic* *asian* *mdhhinc
foreach x of varlist _all {
	rename `x' D_`x'
} 
rename D_COname descnty
gen year=2010
sort year descnty
tempfile tmp
save `tmp', replace
use nimby_estimate, clear
sort year descnty
merge m:1 year descnty using `tmp'
erase `tmp'
drop if _merge==2 //_merge==1 are not 2010, _merge==2 are counties that doesn't generate trash in 2010
drop _merge
replace D_COmdhhinc=D_COmdhhinc*index2000
save nimby_estimate, replace

// 3/3. Import demographics for community of destination
use $data\buffer\demographics_buffer_3mile, clear
keep code year Bpop B0pop Bwhite Bblack Basian Bhispanic Bmdhhinc
local varl "white black asian hispanic"
foreach v in `varl' {
  gen B`v'_p=B`v'*100/Bpop
}
foreach x of varlist _all {
	rename `x' D_`x'
} 
rename D_code Code
rename D_year year
sort year Code
tempfile tmp
save `tmp', replace
use nimby_estimate, clear
sort year Code
merge m:1 year Code using `tmp'
erase `tmp'
drop if _merge==2 //_merge==1 must be 0 if analyzing at other years, _merge==2 are counties that doesn't generate trash in 2010
drop _merge
replace D_Bmdhhinc=D_Bmdhhinc*index2000 if year==2010  // adjust income to 2000 dollars
save nimby_estimate, replace

gen D_Bmdhhinc000=D_Bmdhhinc/1000
sum D_Bpop // check if there is affected community noone lives in
local r = r(min)
if `r'==0 {
replace D_Bmdhhinc=0 if D_Bmdhhinc==.
replace D_Bmdhhinc000=0 if D_Bmdhhinc000==.
replace D_Bblack_p=0 if D_Bblack_p==.
replace D_Bwhite_p=0 if D_Bwhite_p==.
replace D_Bhispanic_p=0 if D_Bhispanic_p==.
replace D_Basian_p=0 if D_Basian_p==.
}
local varv "Operation Access"
foreach v in `varv' {
	tabulate `v', gen (`v')
}
sort yq C J F
save nimby_estimate, replace
******* export to estimate the structural model using MATLAB *********
export delimited using "$matlabin/nimby_estimate", nolabel replace 
 
************ Part of TABLE 2 ************************************************
****** First stage results of the instrumental variable, TABLE 2 ************
lsfereg realprice sumothmktsizes50s isumothmktsizes50s dridis localf, cluster(C) absorb(F T)
test sumothmktsizes50s isumothmktsizes50s

***************************************************************************
************ ROBUSTNESSCHECK ******************************************
**************************************************************************
*** TEST 3.6.2: GENERATE SUBSAMPLES OF DIFFERENT PERIODS: EXIT CHECKS ****
*** RUN ESTIMATION CODE IN MATLAB TO GENERATE TABLES C3 ****
use nimby_estimate, clear
// 15 years
preserve
keep if year>=1995 & year<=2009
export delimited using "$matlabin/nimby_estimate_1995_2009", nolabel replace 
restore 
preserve
keep if year>=2001 & year<=2015
export delimited using "$matlabin/nimby_estimate_2001_2015", nolabel replace 
restore
// 9 years
preserve
keep if year>=2007 & year<=2015
export delimited using "$matlabin/nimby_estimate_2007_2015", nolabel replace 
restore
preserve
keep if year>=1995 & year<=2003
export delimited using "$matlabin/nimby_estimate_1995_2003", nolabel replace 
restore 

********************************************************************
*********** TEST 3.6.3: INDUSTRY STRUCTURE  ********************
use nimby_estimate, clear

replace OperatorEntity=OwnerEntity if OperatorEntity==""
gen firm=OperatorCode
replace firm="Z Public" if OperatorCode=="" & strpos(OperatorEntity,"County")>0 & strpos(OwnerEntity,"County")>0
replace firm="Z Public" if OperatorCode=="" & strpos(OperatorEntity,"City of")>0 & strpos(OwnerEntity,"City of")>0
replace firm="Z Private" if firm==""
replace firm="Z Private" if firm=="NRWS" // small local/1 county company
replace firm="Z Private" if firm=="BURR"
replace firm="Z Private" if firm=="COV"
replace firm="Z Public" if firm=="USDOD"
replace firm="RECO" if firm=="NORCA"
replace firm="Z Private" if firm=="RECO" & OperatorEntity=="Western Placer Recovery Company" // this is a private company Nortech Waste, not Recology, Inc.

egen firmcat=group(firm), label
label define firmcat 8 "Oth private" 9 "Oth public", modify

replace OperatorCode=OwnerCode if OperatorCode==""
replace OperatorCode=OperatorEntity if OperatorCode==""

egen tagtf=tag(yq OperatorCode) if q>0
bys yq: egen nfirm=total(tagtf)
cap gen outlier=dridis>60

************* GRAPH BY COUNTY **********************************************
*number of firms by county: within 60 miles, within, and receive trash
egen tagtfc_w=tag(yq OperatorCode county) if outlier==0
egen tagtfc_wr=tag(yq OperatorCode county) if q>0 & outlier==0
bys yq county: egen nfirmc_w=total(tagtfc_w) // # firms available
bys yq county: egen nfirmc_wr=total(tagtfc_wr) // # firms receive trash
egen tagtlc_w=tag(yq Code county) if outlier==0
egen tagtlc_wr=tag(yq Code county) if q>0 & outlier==0
bys yq county: egen nlandc_w=total(tagtlc_w) // # lands available
bys yq county: egen nlandc_wr=total(tagtlc_wr) // # landfills receive trash
label var nfirmc_w "# firms available"
label var nfirmc_wr "# firms q>0"
label var nlandc_w "# landfills available"
label var nlandc_wr "# landfills q>0"

bys yq county: egen tqc=total(q)
gen shc=tqc/CAq0*100

egen tagtc=tag(yq county)
gsort yq -tagtc shc 
by yq: gen councat_del=_n if  tagtc==1
gen councat_del1=councat_del if year==2015 & quarter==1
bys county: egen councat=mean(councat_del1)
*labmask councat`i' if !missing(councat1) & year==2015& quarter==1, values(county)
labmask councat , values(county)
drop *_del*

replace shc=round(shc,0.1)

********************  FIGURE C2  *******************************************
// # landfills and firms by county
sort yq councat
graph twoway ( bar nlandc_w councat if year==2015 & quarter==1 & councat<=25, fintensity(inten20)) ///
	(bar nfirmc_w councat if year==2015&quarter==1 & councat<=25, fcolor(green) color(green) fintensity(inten50)) ///
	(scatter nlandc_wr councat if year==2015 & quarter==1 & councat<=25, mcolor(navy) msymbol(o) msize(medlarge)) ///
	(scatter nfirmc_wr councat if year==2015 & quarter==1 & councat<=25, mcolor(red) msymbol(x)) ///
	(connected shc councat if year==2015 & quarter==1 & councat<=25, lpattern(solid) lcolor(black*.6) mlabel(shc) mlabposition(9) mlabcolor(black) mcolor(black) msymbol(s) msize(small) yaxis(2)), ///
	xlabel(1(1)25,valuelabel angle(45)) ylabel(#10, angle(45) axis(1)) ///
	 ylabel(#8, angle(45) axis(2))  ///
	legend(off) ///
	xtitle("",placement()) ytitle("# landfills or firms", axis(1))  ytitle("% waste of CA", axis(2)) ///
	graphregion(color(none)) bgcolor(white)
graph export "$result/figureC2_1.pdf", replace

graph twoway ( bar nlandc_w councat if year==2015 & quarter==1 & councat>=26 & councat<=51, fintensity(inten20)) ///
	(bar nfirmc_w councat if year==2015&quarter==1 & councat>=26 & councat<=51, fcolor(green) color(green) fintensity(inten50)) ///
	(scatter nlandc_wr councat if year==2015 & quarter==1 & councat>=26 & councat<=51, mcolor(navy) msymbol(o) msize(medlarge)) ///
	(scatter nfirmc_wr councat if year==2015 & quarter==1 & councat>=26 & councat<=51, mcolor(red) msymbol(x)) ///
	(connected shc councat if year==2015 & quarter==1 & councat>=26 & councat<=51, lpattern(solid) lcolor(black*.6) mlabel(shc) mlabposition(9) mlabcolor(black) mcolor(black) msymbol(s) msize(small) yaxis(2)), ///
	xlabel(26(1)51,valuelabel angle(45)) ylabel(#10, angle(45) axis(1)) ///
	 ylabel(#8, angle(45) axis(2))  ///
	legend(order(1 2 3 4) rows(2) size(*1) symxsize(*0.7)) ///
	xtitle("",placement()) ytitle("# landfills or firms", axis(1))  ytitle("% waste of CA", axis(2)) ///
	graphregion(color(none)) bgcolor(white)
graph export "$result/figureC2_2.pdf", replace

********************* GRAPH BY FIRMS ***********************************
** graph by firms: # counties, # landfills 
*number of counties served by each firm, q>0
egen tagtfic_w=tag(yq firm county) if outlier==0
egen tagtfic_wr=tag(yq firm county) if q>0 & outlier==0

bys yq firm: egen ncounf_w=total(tagtfic_w)
bys yq firm: egen ncounf_wr=total(tagtfic_wr) // # counties receive trash
*number of landfills
egen tagfil=tag(yq firm Code)
bys yq firm: egen nlandf=total(tagfil) // # landfills
bys yq: egen tq=total(q)
bys yq firm: egen tqfi=total(q) // share of CA trash
gen shfi=tqfi/tq*100

egen tagtfi=tag(yq firm)
gsort yq -tagtfi shfi 
by yq: gen firmcat1_del=_n if  tagtfi==1
gen firmcat1_del1=firmcat1_del if year==2015 & quarter==1
bys firm: egen firmcat1=mean(firmcat1_del1)
labmask firmcat1 , values(firm)
drop *_del*
tab firmcat1
label define firmcat1 4 "Oth private" 7 "Oth public", modify

replace shfi=round(shfi,0.1)

label var ncounf_w "# counties within 60 miles"
label var ncounf_wr "# counties send q>0 to"
label var nlandf "# landfills"
label var shfi "% of CA waste"

************************* FIGURE C1-A **********************************
sort yq firmcat1
graph twoway ( bar ncounf_w firmcat1 if year==2015 & quarter==1 , fintensity(inten20)) ///
	(scatter ncounf_wr firmcat1 if year==2015 & quarter==1 , mcolor(navy) msymbol(o) msize(medlarge)) ///
	(scatter nlandf firmcat1 if year==2015 & quarter==1, mcolor(red) msymbol(x)) ///
	(connected shfi firmcat1 if year==2015 & quarter==1, lpattern(solid) lcolor(black*.6) mlabel(shfi) mlabposition(9) mlabcolor(black) mcolor(black) msymbol(s) msize(small) yaxis(2)), ///
	xlabel(1(1)7,valuelabel angle(45)) ///
	ylabel(, axis(2)) ylabel(#10, angle(45) axis(1)) ytitle("% of CA waste", axis(2)) ytitle("# counties or landfills", axis(1)) ///
	legend(rows(4) ring(0) position(10) size(*1) symxsize(*0.7)) ///
	xtitle("",placement())  ///
	graphregion(color(none)) bgcolor(white)
graph export "$result/figureC1a.pdf", replace

*number of landfills and share by county, then average conditional on q>0 county
egen tagtficl_w=tag(yq firm county Code) if outlier==0
egen tagtficl_wr=tag(yq firm county Code) if q>0 & outlier==0

local varv "w wr"
foreach v in `varv' {
bys yq county firm: egen nlandfc_`v'=total(tagtficl_`v')
bys yq firm: egen avgnlandfc_`v'_del=mean(nlandfc_`v') if tagtficl_`v'==1
bys yq firm: egen avgnlandfc_`v'=mean(avgnlandfc_`v'_del)
bys yq firm: egen mednlandfc_`v'_del=pctile(nlandfc_`v') if tagtficl_`v'==1, p(50)
bys yq firm: egen mednlandfc_`v'=mean(mednlandfc_`v'_del)
bys yq firm: egen q1nlandfc_`v'_del=pctile(nlandfc_`v') if tagtficl_`v'==1, p(25)
bys yq firm: egen q1nlandfc_`v'=mean(q1nlandfc_`v'_del)
bys yq firm: egen q3nlandfc_`v'_del=pctile(nlandfc_`v') if tagtficl_`v'==1, p(75)
bys yq firm: egen q3nlandfc_`v'=mean(q3nlandfc_`v'_del)
drop *_del
}

* share by county, then average conditional on q>0 county
cap drop tqc
bys yq county: egen tqc=total(q)
bys yq county firm: egen tqfc=total(q)
gen shfc=tqfc/tqc*100

local varv "w wr"
foreach v in `varv' {
bys yq firm: egen avgshfc_`v'_del=mean(shfc) if tagtfic_`v'==1
bys yq firm: egen avgshfc_`v'=mean(avgshfc_`v'_del)
bys yq firm: egen medshfc_`v'_del=pctile(shfc) if tagtfic_`v'==1, p(50)
bys yq firm: egen medshfc_`v'=mean(medshfc_`v'_del)
bys yq firm: egen q1shfc_`v'_del=pctile(shfc) if tagtfic_`v'==1, p(25)
bys yq firm: egen q1shfc_`v'=mean(q1shfc_`v'_del)
bys yq firm: egen q3shfc_`v'_del=pctile(shfc) if tagtfic_`v'==1, p(75)
bys yq firm: egen q3shfc_`v'=mean(q3shfc_`v'_del)
drop *_del
}
gen firmcat2=firmcat1+0.2
sort firmcat1
************************* FIGURE C1-B **********************************
graph twoway ( scatter mednlandfc_wr firmcat1 if year==2015 & quarter==1, connect(l) mcolor(red*.7) lcolor(red) msymbol(th) ) ///
	(rcap q3nlandfc_wr q1nlandfc_wr firmcat1 if year==2015&quarter==1 , lwidth(thin) lcolor(red) )  ///
	( scatter medshfc_wr firmcat2 if year==2015 & quarter==1, msymbol(th) mcolor(blue*.7) connect(l) lcolor(blue) lpattern(longdash) yaxis(2)) ///
	(rcap q3shfc_wr q1shfc_wr firmcat2 if year==2015&quarter==1 , lwidth(med) lpattern(longdash)  lcolor(blue) yaxis(2) ) , ///
	xlabel(1 "RECO" 2 "ATHEN" 3 "WCNX" 4 "Oth private" 5 "WM" 6 "REPUB" 7 "Oth public", angle(45)) ylabel(#8, axis(2)) ///
	xtitle("",placement()) ytitle("# landfills (with q>0) per county", axis(1))  ///
	ytitle("(%) trash share per county, conditional on q>0", axis(2)) ///
	graphregion(color(none)) bgcolor(white) ///
	legend(order(2 "# landfills per county" 4 "trash share per county (%)" ) rows(2) ring(0) position(10) size(*1) symxsize(*0.7)) 
graph export "$result/figureC1b.pdf", replace
	
*********** EXPORT TO SUBSAMPLES FOR ROBUSTNESSCHECK ****************	
preserve
bys yq county: egen mark=total(inlist(firm,"RECO","REPUB","WCNX","WM"))
replace mark=1 if mark>=1 
keep if mark==1
drop tagtf- mark
sort yq C J F
export delimited using "$matlabin/nimby_estimate_big_firm1", nolabel replace
restore

bys yq county: egen mark=total(inlist(firm,"REPUB","WCNX","WM"))
replace mark=1 if mark>=1 
keep if mark==1
drop tagtf- mark
sort yq C J F
export delimited using "$matlabin/nimby_estimate_big_firm2", nolabel replace


***Using these delimited files as inputs to estimate costs and counterfactual...
//... waste flows in MATLAB
erase index2000.dta
erase nimby_estimate.dta
erase nimby_main.dta























